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Abstract: Binary sequences with low autocorrela- 
tions are important in communication engineering and in 
statistical mechanics as groundstates of the Bernasconi- 
modcl. Computer searches are the main tool to construct 
such sequences. Due to the exponential size 0(2 ) of 
the configuration space, exhaustive searches are limited to 
short sequences. We discuss an exhaustive search algo- 
rithm with run time characteristic 0(1.85^) and apply it 
to compile a table of exact groundstates of the Bernasconi- 
model up to N — 48. The data suggests F > 9 for the 
optimal merit factor in the limit N — > oo. 

1 Introduction 



Binary sequences S ■■ 
autocorrelations 



{si = ±1, . . . , sn} with low off-peak 



C k (S) 



N-k 



(1) 



i=l 



have applications in many communication engineering 
problems [1]. One exciting example has been their use in 
high precision interplanetary radar measurements to check 
out space-time-curvature [2]. 

Physicists prefer to consider binary sequences as one 
dimensional systems of Ising-spins. In this context, low 
autocorrelation binary sequences appear as minima of the 
energy 



JV-l 



e{s) = C K S )- 



(2) 



fe=i 



This is the Bernasconi- model [3] . It has long-range 4-spin 
interactions and is completely deterministic, i.e. there is 
no explicit or quenched disorder like in spin-glasses. Ne- 
vertheless the ground states are highly disordered - quasi 
by definition. This self-induced disorder resembles very 
much the situation in real glasses. In fact, the Bernasconi- 
modcl exhibits features of a glass transition like a jump in 
the specific heat [3] and slow dynamics and aging [4] . 

A clever variation of the replica method allows an ana- 
lytical treatment of the Bernasconi-model in the high- 



temperature regime [5, 6]. For the low-temperature re- 
gime, analytical results are rare - especially the ground 
states are not known. With periodic boundary conditions, 
i.e. with 

JV 

Cfc = ^ s i s (i+k-l)( mod N)+l> (3) 
i=l 

instead of eq. 1, the construction of ground states is pos- 
sible for special values of N. Example: For N = An + 3 
prime, the modified Legendre-scqucncc 

J j|(iV-l) ft l<i<N ( .s 

S i _ 1 ±1 A - AT [ 4:1 



N 



yields C% = 1, the minimum possible value for odd N. 
Other ground states can be constructed from linear shift 
register sequences based on primitive polynomials over Ga- 
lois fields. This construction requires N — 2 P — 1 with p 
prime. See [1, 6] for details. 

For the model with open boundary conditions (eq. 1) 
no construction of ground states is known, not even for 
special values of N. The Legendre-sequences are far from 
the true ground states [7]. The only exact results have 
been provided by exhaustive enumerations, limited how- 
ever by the exponential complexity of the problem to sys- 
tems smaller than N — 32 [8]. Partial enumerations al- 
low larger values of N but cannot guarantee to yield true 
ground states. Promising candidates for partial enumer- 
ations are the skew-symmetric sequences of odd length 
N = 2n — 1. These sequences satisfy 



Sn+l = (-1) 



1 = 1,. 



1 



(5) 



from which it follows that all Ck with k odd vanish. The 
restriction to skew-symmetric sequences reduces the effec- 
tive size of the problem by a factor of 2, but the true 
ground states are not skew-symmetric for several values of 
N, as we will see below. 

Finding the ground states of the Bernasconi-model has 
turned out to be a hard mathematical problem. Golay 
[8, 3] has conjectured that the maximal merit factor 
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should obey F < 12.32 for N » 1. However, heuristic 
searches among skew-symmetric sequences up to TV = 199 
suggest F « 6 for long sequences [9], a value consis- 
tent with results from simulated annealing [3]. This large 
discrepancy indicates that the groundstates, i.e. the se- 
quences with high merit factors 6 < F < 12, must be 
extremely isolated energy minima in configuration space. 
Stochastic search procedures including simulated anneal- 
ing are not well suited to find these "golf-holes". Ex- 
haustive search seems to be the only approach at least for 
small systems. The complete configuration space has been 
searched up to N = 32 [8], the skew-symmetric subspace 
up to N = 71 [10, 11]. Fifty days of CPU-time on a spe- 
cial purpose computer have been used for an exhaustive 
search for binary sequences up to N = 40 that minimize 
max fe \ C k \ [12]. 

In this contribution we discuss a fast algorithm for the 
exhaustive enumeration. It is fast enough to yield exact 
groundstates of the Bernasconi-model up to N — 48 and 
can easily be modified for partial enumerations. The data 
is used to estimate the optimal merit factor in the large 
N limit. 



2 The Algorithm 

Any algorithm that performs an exhaustive search for the 
ground state of the Bernasconi-model has to cope with 
the enormous size (2^) of the configuration space. This 
exponential complexity limits the accessible values of N 
very soon and calls for methods to restrict the search to 
smaller subspaces without missing the true ground state. 
Symmetries are an abvious device to cut out portions of 
the configuration space. We will see below, that the use 
of symmetries can reduce the size of the search space by a 
factor of about 1/8. A method borrowed from combinato- 
rial optimization - branch and bound - will prove useful to 
reduce the complexity from 0(2 N ) to 0(b N ) with b < 2. 
We will further see, that the enumeration problem is suited 
almost perfectly for parallelization. 

2.1 Symmetries 

The correlations Ck (eq. 1) are unchanged when the se- 
quence is complemented or reversed. When alternate ele- 
ments of the sequence are complemented, the even-indexed 
correlations are not affected, the odd-indexed correlations 
only change sign. Hence, with the exception of a small 
number of symmetric sequences, the 2 N sequences will 
come in classes of eight which are equivalent. The total 
number of nonequivalent sequences is slightly larger than 
2 N - 3 . 

The to left- and to rightmost elements of the sequence 
can be used to parameterize the symmetry-classes. For 



to = 3 and N odd, this gives 12 classes: 
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For N even there are 10 classes. In general the number c 
of symmetry-classes that can be distingiushed by to left- 
and to right-border elements reads 

c(m) = 2 2m - 3 + 2 m - 2+{N mod 2) (7) 

and the number of nonequivalent configurations reduces 
to a fraction 

< m ) = 1 1 m 

22m g "'" 2m+2-(Af mod 2) ' ^ ' 

The optimal value | is approached with increasing to. 
2.2 Branch and Bound 

Branch and bound methods are commonly used in combi- 
natorial optimization [13] and (less frequently) in statisti- 
cal mechanics [14, 15]. They solve a discrete optimization 
problem by breaking up its feasible set into successively 
smaller subsets {branch), calculating bounds on the ob- 
jective function value over each subset, and using them to 
discard certain subsets from further consideration {bound). 
The procedure ends when each subset has either produced 
a feasible solution, or has been shown to contain no better 
solution than the one already in hand. The best solution 
found during this procedure is a global optimum. 

The idea is of course to discard many subsets as early as 
possible during the branching process, i.e. to discard most 
of the feasible solutions before actually evaluating them. 
The success of this approach depends on the branching 
rule and very much on the quality of the bounds, but 
it can be quite dramatic. Numerical investigations have 
shown, for example, that the n-city Traveling Salesman 
Problem can be solved exactly in time 0{n a ) with a < 3 
using branch and bound methods [13]! This is no contra- 
diction to the exponential complexity of the TSP since the 
latter is the guaranteed, i.e. worst case complexity, while 
the former refers to the typical case, averaged over many 
instances of the TSP. 

In accordance with our symmetry-classes, we specify a 
set of feasible solutions by fixing the to left- and to right- 
most elements of the sequence. The N — 2to center ele- 
ments are not specified, i.e. the set contains 2 JV ~ 2m fea- 
sible solutions. Given a feasible set specified by the to 
border elements, 4 smaller sets are created by fixing the 
elements s m+ i and s^^ m to ±1. This is the branching 
rule. It is applied recursively until all elements have been 
fixed. The energy of the resulting sequence is compared 
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Algorithm 1 Procedure search (S, m) - search for the 
minimum energy configuration S opt within the subset (S, 
to) of all configurations. 
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n <— N — 2m: {number of free elements Si} 
if n > 1 then {> 2 seq. in subset} 
if E b (S,m) > E(S opt ) then {bound} 

return; 
else {branch} 

search (S, m + 1); 

s m +i < Sm+i; search (5, m + 1); 

SiV-m < SAT_ m ; search (5, to + 1); 

Sm+i < «m+i; search (S, to + 1); 

end if 

else if n = 1 then {2 seq. in subset} 
if £(S) < E(S opt ) then 

Sopt < S 1 ; 
end if 

if E(S) < E(S opt ) then 

$opt * S ; 
end if 
else {1 seq. in subset} 
if E(S) < E(S opt ) then 

Sopt < S 1 ; 
end if 
end if 



to the minimum energy found so far. If it is lower, the 
sequence is kept as the potential ground state. After all 
c(m)2 N ~ 2m sequences have been considered, the potential 
ground state has turned into a true one. 

Lower bounds are usually obtained by replacing the 
original problem over a given subset with an easier (re- 
laxed) problem such that the solution value of the latter 
bounds that of the former. A good relaxation is one that: 
(i) is easy and fast to solve; and (ii) yields strong lower 
bounds. Most often these are conflicting goals. 

A relaxation of the LABS problem is given by the prob- 
lem to adjust the free elements (i.e. the center elements 
s m+ i, . . . , sjv_ m ) to minimize all values C| separately, i.e. 
we replace the original problem 



'N-l 

Emin = min | ^ C l 



subset 



, fc=l 



with the relaxed version 

SU DSCt 



fe=i 



(9) 



(10) 



Unfortunately E^ in is still not easy to calculate, but we 
can proceed our relaxation by providing a lower bound 
Eb < -E^in- For that purpose we choose an arbitrary se- 
quence from the subset with correlations Ck- Comple- 
menting a free element Si — Si can decrement \Ck\ at 



most by 2. Let fk denote the number of terms SiS i+ k in 
Ck that contain at least one free element. This leads to 

N-k 

E b = ^min{& fc ,(|C fe |-2/ fe ) 2 } (H) 
fe=i 

< £min < ^min 

where b k = (N — k) mod 2 6 {0, 1} is the minimum value 
\Ck\ can attain. The fk are given by 



fk={ 2(7V- 
N- 





- TO — k) 

2m 



k>N-m 

N/2 < k < N - to 

k < N/2 



(12) 



i.e. the long range correlations are not affected by our 
relaxation. E b is not the strongest bound to E m i n , but 
its calculation is very fast. 

Now we have gathered all ingredients to formulate the 
branch and bound procedure search (algorithm 1). This 
procedure is called with two parameters specifying the 
subset to search: a binary sequence S — {si, . . . , sn} and 
an integer m. The subset consists of all 2 N ~ 2m sequences 
that can be generated from S by varying the N — 2m 
center elements. 5 op t is a global variable that holds the 
sequence with the minimum energy found so far. On entry, 
the size of the subset is checked: If it contains more than 
2 sequences, branch and bound (lines 3-10) is applied. 
Otherwise the sequences in the subset are evaluated (lines 
11-23). 

The procedure search is called from a driving proce- 
dure with c(m ) subsets, each representing a different 
symmetry-class. In practice, we used mo = 6 with 528 
(N even) resp. 544 (N odd) symmetry-classes. 

Figure 1 CPU time for exhaustive search algorithm vs. 
N. Times are measured on a Sun UltraSparc I 170 work- 
station. 
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To measure the impact of branch and bound, we started 
two runs on the same machine: one "straight" enumer- 
ation (omitting lines 3-5) and the other with activated 
bound-mechanism. Figure 1 displays the effect of branch 
and bound on the CPU time. The straight enumeration 
shows the expected 0(2^) behavior. Branch and bound 
reduces the complexity to 0(1. 85^). Albeit this is still ex- 
ponential, the gain in speed is worth the little effort. The 
branch and bound enumeration for N = took about 2 
days on a Sun UltraSparc 1 170 workstation. This com- 
pares well with the extrapolated 68 days for the straight 
enumeration! 

2.3 Parallelization 

The different symmetry classes can be searched indepen- 
dently. Hence the straight enumeration is perfectly paral- 
lclizable into c(m) threads of control. Branch and bound 
complicates the situation: Whenever a better sequence is 
found by one thread, it should be communicated immedi- 
ately to all other threads to ensure that always the best 
E(S opt ) is used in the bounding test (line 3). But E(S opt ) 
is accessed very frequently, so the necessary synchroniza- 
tion would spoil the parallelization. Giving each thread 
its own local copy of S op t preserves perfect parallelization 
but abandons most of the benefits of branch and bound! 



Figure 2 Speed-up factor of the branch and bound algo- 
rithm on a symmetric multiprocessor platform using 2, 3 
or 4 CPUs. 



new work assignment. This process repeats until all sym- 
metry classes have been considered. 

The access to the workpile has to be protected with a 
mutual exclusion lock, allowing only one thread at a time 
to read or modify data from the workpile. If each worker 
thread uses its own local copy of S op t, this is the only syn- 
chronization needed. To propagate the best S'opt as fast as 
possible among the workers, it is stored in the workpile. 
A worker that requests a new work assignment compares 
its own local S op t with the global one and updates the 
one with the higher energy under the protection of the 
lock. This method limits the use of a suboptimal S'opt to 
the search within n — 1 symmetry classes, where n is the 
number of worker threads. The delay in propagation of 
the optimal S op t is minimized by choosing c(m) ^> n. In 
this case, the workpile paradigm has the additional ad- 
vantage of evenly distributing the load among all worker 
threads. Due to branch and bound, the enumerations in 
some symmetry classes may take considerably less time 
than in others. A worker that encounters these "easy" 
classes simply gets more classes to search. 

On a 4 processor Sun SPARCstation 20, the number 
of worker threads (< 4) is much smaller than c(m) for 
m = 6, so the workpile paradigm should yield almost per- 
fect parallelization. Figure 2 shows that this is indeed the 
case. The low speed-up factors for small N are due to 
the relative costs of thread-generation and synchroniza- 
tion compared to the actual enumeration. 

3 Results 
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A solution to this dilemma is provided by the workpile 
paradigm [16]: The symmetry classes to be searched for 
are put on a central workpile and a number of worker 
threads are launched together. Each worker thread re- 
quests an assignment of work from the workpile (i.e. a 
symmetry class), performs the search, and then asks for a 



Using the multithreaded branch and bound algorithm and 
313 hours of CPU time on a 4-processor Sun SPARCsta- 
tion 20, the ground states of the Bernasconi-model have 
been found up to N = 48 (table 1). The enumeration for 
N = 32 (the previous peak value) took only 80 seconds, 
N — 39 was done in 1 hour. It is remarkable that from 
the 22 optimal skew-symmetric sequences in the range 
5 < N < 47 [10] 7 have energies well above the true 
ground-state energie - a third. This should be kept in 
mind if one uses skew-symmetric sequences to estimate 
the groundstate energy in the limit N — > oo. 

Figure 3 shows the groundstate energies E vs. N. In 
contrast to the model with periodic boundary conditions 
there arc no visible regular patterns for special values of 
N [6]. The energies seem to follow E oc N 2 for all values 
of N. A quadratic fit yields 



F 



_/V 2 
hm 

N->oc 2E 



9.3 



and leads us to the tentative conclusion that 

2E 



F= lim 



> 9. 



(13) 



(14) 



This estimate is in agreement with Golay's conjecture F < 
12.32 and has to be compared to the value F « 6.0 found 
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Table 1 Ground states of the Bernasconi-model for 3 < 
N < 48. Sequences are written in run- length notation: 
Each figure indicates the number of consecutive equal 
signed elements. 
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E ■ 
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21 
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211 


5 


2 


311 


6 
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1113 


7 


3 


1123 


8 


8 


12113 


9 


12 


42111 


10 


13 


22114 


11 


5 


112133 


12 


10 


1221114 


13 


6 


5221111 


14 


19 


2221115 


15 


15 


52221111 


16 


24 


225111121 


17 


32 


252211121 


18 


25 


441112221 


19 


29 


4111142212 


20 


26 


5113112321 


21 


26 


27221111121 


22 


39 


51221111233 


23 


47 


212121111632 


24 


36 


2236111112121 


25 


36 


337111121221 


26 


45 


21212111116322 


27 


37 


34313131211211 


28 


50 


34313131211212 


29 


62 


212112131313431 


30 


59 


551212111113231 


31 


67 


7332212211112111 


32 


64 


71112111133221221 


33 


64 


742112111111122221 


34 


65 


842112111111122221 


35 


73 


7122122111121111332 


36 


82 


3632311131212111211 


37 


86 


844211211111122221 


38 


87 


8442112111111122221 


39 


99 


82121121234321111111 


40 


108 


44412112131121313131 


41 


108 


343111111222281211211 


42 


101 


313131341343112112112 


43 


109 


1132432111117212112213 


44 


122 


525313113111222111211121 


45 


118 


82121121231234321111111 


46 


131 


823431231211212211111111 


47 


135 


923431231211212211111111 


48 


140 


3111111832143212221121121 



Figure 3 Groundstate energy of the Bernasconi- 
Hamiltonian vs. N. 
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by heuristic searches for long skew-symmetric sequences [9] 
and by simulated annealing [3]. This indicates once more 
that heuristic and probabilistic methods fail to find the 
groundstates of the Bernasconi-model. Every algorithm of 
such a kind should be jugded by the percentage of values 
it finds from table 1. 
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